disp('fmdlf');

% disp(fn_solve(33.0));
% 
% Num = 100;
% h0 = 1.0;
% s = 1.0;
% W = h0 + 2*s;
% vy = fn_solve_vy(Num, h0, s);
% 
% f = @(x) log(1+exp(-x));
% 
% for n = 1:length(vy)
%     
%    tmp = vy(n) + f(vy(n));
%    for n1 = 1:(n-1)
%        tmp = tmp + 2*(n - n1)*f(vy(n1));
%    end
%    
%    tmp = tmp - n*W;
%    fprintf('n = %d, diss = %e \n', n, tmp);
%     
% end

% disp(vy');

Num = 400;
h0 = 0.1;
vs0 = linspace(0.686, 0.9, 400);
vs = zeros(1, length(vs0));
vexps = zeros(1, length(vs0));
for i = 1:length(vs0)
    s0 = vs0(i);
    s = fn_s_to_s(Num, h0, s0);
    vs(i) = s;
end

figure(8); plot(vs0, vs, '-ro');
figure(9); plot(vs0, vs0 - vs, '-bo');

[s_left, s_right] = fn_total_solve_lr(Num, h0);
[vy, s_mid] = fn_total_solve(Num, h0);
fprintf('s_mid_judge = %e \n', fn_judge(Num, h0, s_mid));

v_diss = zeros(1, length(vy));

f = @(x) log(1+exp(-x));
for n = 1:length(vy)
    
   tmp = vy(n) + f(vy(n));
   for n1 = 1:(n-1)
       tmp = tmp + 2*(n - n1)*f(vy(n1));
   end
   
   tmp = tmp - n*h0 - 2*n*sum(log(1+exp(-vy)));
%    fprintf('n = %d, diss = %e \n', n, tmp);
   tmp_real = real(tmp);
   tmp_imag = imag(tmp);
   tmp_imag = mod(tmp_imag+pi, 2*pi) - pi;
   tmp = tmp_real + 1j*tmp_imag;
   v_diss(n) = tmp;
    
end

fprintf('max loss = %e \n', max(abs(v_diss)));

% close(figure(6));

disp('fmdlf');

Num = 800;
vy2 = fn_total_solve(Num, h0);

close(figure(99));
figure(99);
plot(1:length(vy), vy, '-r*');
hold on;
plot(1:length(vy2), vy2, '-bo');

disp('fmdkx');

% Num = 1000;
% h0 = 0.;
% 
% [s_left, s_right] = fn_total_solve_lr(Num, h0);
% vy = fn_total_solve(Num, h0);
% 
% v_diss = zeros(1, length(vy));
% 
% 
% for n = 1:length(vy)
%     
%    tmp = vy(n) + f(vy(n));
%    for n1 = 1:(n-1)
%        tmp = tmp + 2*f(vy(n1));
%    end
%    
%    tmp = tmp - n*h0 - 2*n*sum(log(1+exp(-vy)));
% %    fprintf('n = %d, diss = %e \n', n, tmp);
%    v_diss(n) = tmp;
%     
% end
% 
% fprintf('max diss = %e \n', max(abs(v_diss)));
% figure(1); hold on;
% plot(1:length(vy), vy, '-o');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function res = fn_find_singularity(Num, h0)

s_left = 0.1;
while ~fn_if_singu(Num, h0, s_left)
    s_left = s_left - 0.1;
end

s_right = 0.9;
while fn_if_singu(Num, h0, s_right)
    s_right = s_right + 0.1;
end

while 1
    
    if s_right - s_left < 1e-12
        break;
    end
    
    s_mid = 1/2*(s_left + s_right);
    if fn_if_singu(Num, h0, s_mid)
        s_left = s_mid;
    else
        s_right = s_mid;
    end

end

res = s_right;

end




function if_singu = fn_if_singu(Num, h0, s0)

vy = fn_solve_vy(Num, h0, s0);
tmp = sum(abs(imag(vy)));
if tmp>0.1
    if_singu = (1>0);
else
    if_singu = (1<0);
end

end



function [vy, s_mid] = fn_total_solve(Num, h0)

[s_left, s_right] = fn_total_solve_lr(Num, h0);

if fn_judge(Num, h0, s_left) == 0
    vy = fn_solve_vy(Num, h0, s_left);
    return;
end

if fn_judge(Num, h0, s_right) == 0
    vy = fn_solve_vy(Num, h0, s_right);
    return;
end

while 1

    s_mid = 1/2*(s_left + s_right);
    
    if abs(s_left - s_right) < 1e-14
        break;
    end
    
    if fn_judge(Num, h0, s_mid) == 0
        vy = fn_solve_vy(Num, h0, s_mid);
        return;
    elseif fn_judge(Num, h0, s_mid)>0
        s_right = s_mid;
    else
        s_left = s_mid;
    end
    
end

vy = fn_solve_vy(Num, h0, s_mid);

end



function [s_left, s_right] = fn_total_solve_lr(Num, h0)

s_left = fn_find_singularity(Num, h0);

s_right = s_left;

while 1
    s_right = s_right + 0.1;
    if fn_judge(Num, h0, s_right) > 0
        break;
    end
end

s_left = s_right - 0.1;

fprintf('s_left = %f, s_right = %f \n', s_left, s_right);
fprintf('judge_left = %f, judge_right = %f, \n', ...
    fn_judge(Num, h0, s_left), fn_judge(Num, h0, s_right));

end



function res = fn_judge(Num, h0, s)

res = s - real(fn_s_to_s(Num, h0, s));

end


function exps = fn_s_to_exps(Num, h0, s0)

vy = fn_solve_vy(Num, h0, s0);
exps = fn_exps(vy);

end


function [s, vy] = fn_s_to_s(Num, h0, s0)

vy = fn_solve_vy(Num, h0, s0);
s = fn_s(vy);

end


function exps = fn_exps(vy)

exps = prod(1 + exp(-vy));

exps = real(exps);

end


function s = fn_s(vy)

s = sum(log(1 + exp(-vy)));

end


function vy = fn_solve_vy(Num, h0, s)

W = h0 + 2*s;
vy = zeros(1, Num);

for n = 1:Num
    
    tv = vy(1:n-1);
    RHS = n*W - 2*sum(((n-1):-1:1).*log(1 + exp(-tv)));
    vy(n) = fn_solve(RHS);
    
end



end



function res = fn_solve(f)

if f < 33
    res = log(exp(f) - 1);
else
    res = f;
end


end